Series Temporales

Introducción

Los datos de este apartado de la práctica son datos de la demanda eléctrica b.c en España. Es decir, a partir de toda la energía que se ha genera, se ha hace el balance de aquella parte que se vende, compra o utiliza para bombero, y finalmente queda aquella que queda disponible para proveer a la población de España.

Ésta ha sido obtenido a través de la API de Red Eléctrica (REE). Se ha adjuntado el fichero con el que se han obtenido los datos. Éstos son diarios y van desde enero 2018 hasta febrero 2020 para que no coincida el periodo que se quiere predecir con el COVID19 ya que dificultaría la práctica.

Lectura de datos

setwd("~/Escritorio/MASTER/CUATRI3/MACHINE_LEARNING2/PRACTICA/ML2/Series Temporales")
df_ree <- read.csv("~/Escritorio/MASTER/CUATRI3/MACHINE_LEARNING2/PRACTICA/ML2/Series Temporales/ree.csv")

df_ree$id <- NULL
df_ree$date = substr(df_ree$datetime,1,10)

# Fucniones
mape <- function(actual,pred){
  mape <- mean(abs((actual - pred)/actual))*100
  return (mape)}

# FUNCIÓN MY RESIDUALS (para comparar bien ACF y)
my_tsresiduals <- function(data, ...){
  if(!fabletools::is_mable(data)){
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  
  data <- stats::residuals(data)
  if(n_keys(data) > 1){
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  
  gg_tsdisplay(data, !!sym(".resid"), plot_type = "partial", ...)
}

#' @export
`+.gg_tsensemble` <- function(e1, e2){
  e1[[1]] <- e1[[1]] + e2
  e1
}

#' @export
print.gg_tsensemble <- function(x, ...){
  print(x[[1]], vp = grid::viewport(layout.pos.row = c(1, 1), layout.pos.col = c(1, 2)))
  print(x[[2]], vp = grid::viewport(layout.pos.row = 2, layout.pos.col = 1))
  print(x[[3]], vp = grid::viewport(layout.pos.row = 2, layout.pos.col = 2))
}

Convertimos la fecha a un formato adecuado.

df_demanda <- df_ree %>%
  filter(subcategory == "Demanda en b.c.")

df_demanda$datetime <- as.POSIXct(df_demanda$datetime)
df_demanda$date <- as.Date(as.character(df_demanda$datetime)) #, format="%d/%m/%Y", origin = "2018-01-01"

Dividimos Train - Test:

• Train: enero 2018- enero 2020
• Test: 1-10 de febrero 2020

df_demanda_train = df_demanda[df_demanda$date < '2020-02-01',]

df_demanda_test = df_demanda[(df_demanda$date > '2020-01-31') & (df_demanda$date <= '2020-02-10'),] 

df_demanda = df_demanda[(df_demanda$date < '2020-02-10'),]

Análisis previo series temporales

Convertimos el dataframe a un objeto tsibble:

df_demanda_train <- df_demanda_train %>%
  as_tsibble(index = date) 

df_demanda_test <- df_demanda_test %>%
  as_tsibble(index = date) 

df_demanda <- df_demanda %>%
  as_tsibble(index = date)

Representamos la serie temporal con la que se trabajará. A priori parece que sí hay cambios sistemáticos en media (no hay tendencia a largo plazo) pero no lo parece en varianza (las líneas de la envolvente en los límites superior e inferior de la serie varían).

df_demanda %>%
  autoplot(energy) +
  labs(title = "Evolución de la demanda eléctrica en España") 

Observamos la descomposición aditiva de la serie. Observamos la serie descompuesta en la componente de tendencia, estacional y el residuo tras eliminar las componentes estacional y de tendencia.

df_demanda %>%
  model(classical_decomposition(energy, type = "additive")) %>%
  components() %>%
  autoplot() 

Representamos la serie diferenciada.

Se observa que tras la diferenciación la serie pierde la tendencia, parece estacionaria (sin cambios sistemáticos en media ni varianza).

df_demanda_train %>% autoplot( difference(energy)) + ggtitle("Diferenciación de la serie de la demanda energética  ")

Para determinar de manera objetiva si necesario diferenciar se utiliza una prueba de raíz unitaria. En este contraste, la hipótesis nula es que los datos son estacionarios, y buscamos pruebas de que la hipótesis nula es falsa.

En consecuencia, como es el caso, pequeños p-valores (por ejemplo, menos de 0,05) sugieren que es necesario diferenciar. El p-valor es menor que 0.05, lo que indica que la hipótesis nula es rechazada:

df_demanda_train %>%
  features(energy, unitroot_kpss)

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1     0.467      0.0491

Vemos que la serie diferenciada obtiene un p-valor mayor que 0.05:

df_demanda_train %>%
  mutate(diff_close = difference(energy)) %>%
  features(diff_close, unitroot_kpss)

## # A tibble: 1 x 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0468         0.1

Autocorrelaciones

Como hemos visto, se trata de una serie no estacionaria, por lo que habrá que aplicar un modelo ARIMA o SARIMA. Este último permite incluir efectos de cambios repetitivos en periodos concretos de la serie, por lo que puede sernos útil.

Para empezar a definir los parámetros p,d,q utilizaremos las gráficas de correlogramas ACF y PACF.

Autocorrelación parcial

Vemos que sigue un mismo patrón sinusoidal (desproporcional) y que va progresivamente disminuyendo. Vemos que tiene un pico alto en el primer retardo y que posteriormente se repiten picos cada 7 retardos. Esto implicaría que tiene una autocorrelación importante con los datos del día anterior y los de hace una semana, éstos últimos indican una componente estacional.

Esto tiene sentido ya que la demanda eléctrica depende mucho de los hábitos de consumo de los particulares, así como de las producciones de las empresas (entre semana y fin de semana).

El pico significativo en el retardo 7 del ACF sugiere un componente MA(1) estacional. Los dos picos significativos en el retardo 1 y 2 del ACF sugiere por lo menos una componente MA(2) no estacional.

df_demanda_train %>% ACF(energy) %>% autoplot()

Autocorrelación parcial

En este caso los valores están más acotados. Sin embargo, se continúan observando picos estacionales a los 7 días. Los dos picos significativos en el retardo 1 y 2, una componente MA(1) o MA(2) no estacional. El pico en el retardo 7 sugiere una componente MA(1) estacional.

df_demanda_train %>% PACF(energy) %>% autoplot()

Como hemos comprobado en ambos gráficos hay una componente estacional importante, por lo que se va a proceder a probar directamente un modelo SARIMA.

SARIMA

Ajuste 1

Como hemos comentado, empezaremos ajustando un modelo que represente las componentes correladas con un número pequeño de valores anteriores (1 o 2). Además, reflejaremos la componente estacional.

fit <- df_demanda_train %>%
  model(arima = ARIMA(energy ~ pdq(0,0,2) + PDQ(0,0,1))) #pdq(0,0,2) + PDQ(0,0,2)

Como se observa, el modelo ya detecta una estacionalidad de 7 días:

report(fit)

## Series: energy 
## Model: ARIMA(0,0,2)(0,0,1)[7] w/ mean 
## 
## Coefficients:
##          ma1     ma2    sma1    constant
##       0.8185  0.2267  0.6022  731630.136
## s.e.  0.0340  0.0354  0.0242    4596.651
## 
## sigma^2 estimated as 1.518e+09:  log likelihood=-9123.85
## AIC=18257.71   AICc=18257.79   BIC=18280.88

El modelo tiene que ser estacionario invertible. Representamos las raíces inversas para comprobar que todas ellas están dentro del círculo unitario.

gg_arma(fit)

Si mostramos los residuos, vemos que aún no es un ruido blanco, aparecen picos que habría que evitar. Si observamos el gráfico ACF vemos que aún hay picos en retardos 14, 21, 28.. que hay que ajustar. Si observamos el gráfico PACF vemos que hay retardos en los mismos retardos con picos prácticamente.

fit %>% my_tsresiduals()

Ajuste 2

Probamos ajustando valores para modelar las componentes estacionales que hemos visto en los gráficos de autocorrelaciones de los residuos. Para ello añadimos D=1 para intentar recoger esa diferenciación semanal.

fit2 <- df_demanda_train %>%
  model(arima = ARIMA(energy ~  pdq(0,0,2) + PDQ(0,1,1))) 

report(fit2)

## Series: energy 
## Model: ARIMA(0,0,2)(0,1,1)[7] 
## 
## Coefficients:
##          ma1     ma2     sma1
##       0.9651  0.3648  -0.7170
## s.e.  0.0342  0.0293   0.0421
## 
## sigma^2 estimated as 845137331:  log likelihood=-8820.64
## AIC=17649.27   AICc=17649.33   BIC=17667.77

Comprobamos que se mantienen las raíces dentro del círculo:

gg_arma(fit2)

Comprobamos los residuos. Vemos que muchos obtienen valores más cercanos a 0 pero quedan muchos aún con valores más altos, pueden corresponder a picos puntuales que quizás se mejorarían añadiendo una variable exógena con las festividades al año.

Sin embargo, en cuanto a los gráficos de las autocorrelaciones vemos que han desaparecido los picos a los 7 días pero han aparecido otros en los retardos 2,3.

fit2 %>% my_tsresiduals()

Ajuste 3

Intentamos reducir las correlaciones en los primeros retardos añadiendo d=1 en la parte no estacional.

fit3 <- df_demanda_train %>%
  model(arima = ARIMA(energy ~  pdq(0,1,2) + PDQ(0,1,1))) 

report(fit3)

## Series: energy 
## Model: ARIMA(0,1,2)(0,1,1)[7] 
## 
## Coefficients:
##           ma1      ma2     sma1
##       -0.0442  -0.3114  -0.9776
## s.e.   0.0362   0.0374   0.0143
## 
## sigma^2 estimated as 608744175:  log likelihood=-8693.38
## AIC=17394.77   AICc=17394.82   BIC=17413.26

Sin embargo con esta configuración las raíces aparecen en los límites del círculo, por lo que no sería invertible y no es recomendable.

gg_arma(fit3)

Aún así se muestran los residuos obtenidos.

fit3 %>% my_tsresiduals()

Se ha intentado modificar más los parámetros pero no se ha logrado obtener una configuración que redujera más los picos de los correlogramas cumpliendo que sea invertible y estacionario. Por ello finalmente se ha procedido a ajustar un modelo SARIMA automáticamente.

Ajuste 4

Se aplica el modelo SARIMA automático.

fit4 <- df_demanda_train %>%
  model(arima = ARIMA(energy))

report(fit4)

## Series: energy 
## Model: ARIMA(2,0,3)(1,1,0)[7] 
## 
## Coefficients:
##          ar1      ar2      ma1      ma2     ma3     sar1
##       1.8261  -0.8760  -0.9543  -0.2272  0.3341  -0.3956
## s.e.  0.0358   0.0333   0.0529   0.0466  0.0484   0.0363
## 
## sigma^2 estimated as 807256041:  log likelihood=-8800.02
## AIC=17614.05   AICc=17614.2   BIC=17646.42

gg_arma(fit4)

Si visualizamos los residuos, tienen un comportamiento parecido al modelo 2 ajustado. Sin embargo, como diferencia, vemos que los picos que hay todavía en el correlograma están en el retardo 14, y posteriormente cada 7.

fit4 %>% my_tsresiduals()

Entre este modelo y el modelo ajuste 2, se ha decidido realizar las predicciones con este ya que hay un número menor de picos y son más tardíos.

Forecast SARIMA

Realizamos la predicción para los 10 días en febrero.

forecast =  forecast(fit4, h=10) 

forecast %>% autoplot(df_demanda%>% filter(date >= "2020-01-01"))

Evaluación SARIMA

Finalmente calculamos el error de la predicción:

accuracy(forecast$energy, df_demanda_test$energy)

##              ME     RMSE   MAE     MPE   MAPE
## Test set -31089 37356.54 31089 -4.2756 4.2756

Vemos que el error es aceptable para el horizonte contemplado.

Modelo Dummy

Para poder evaluzar el modelo y compararlo, se calcula un modelo dummy Naive estacional. Siempre aporta valor tener un punto de referencia. Para ello entrenamos en un mismo objeto el SARIMA obtenido antes y el Naive.

fit_dummy <- df_demanda_train %>%
  model(
    arima = ARIMA(energy ~  pdq(2,0,3) +  PDQ(1,1,0)),
    snaive = SNAIVE(energy)) 

Realizamos la predicción:

fc <- fit_dummy %>% forecast(h = 10)

Representamos los resultados:

fc %>%
  filter(date >= "2020-01-01") %>% autoplot(df_demanda %>% filter(date >= "2020-01-01"))

Calculamos finalmente los errores:

accuracy(fc, df_demanda_test)

## # A tibble: 2 x 9
##   .model .type      ME   RMSE    MAE   MPE  MAPE  MASE  ACF1
##   <chr>  <chr>   <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 arima  Test  -31089. 37357. 31089. -4.28  4.28   NaN 0.296
## 2 snaive Test  -39454. 46292. 40951. -5.59  5.79   NaN 0.504

Observamos que el MAPE del modelo SARIMA es un 1,3% menor que el modelo Naive, por lo que consideramos que es mejor el SARIMA.

Conclusiones

Como conclusiones queríamos comentar que se han probado muchas configuraciones del modelo SARIMA, muchas no eran estables o no cumplían las condiciones de ser invertibles y estacionarias, tampoco acababan de cumplir todas las condiciones para unos residuos de ruido blanco, aunque no están muy alejadas.

Otra opción que se podría analizar es realizar un modelo de series temporales que prediga sólo la demanda del siguiente día y se vaya actualizando con los datos reales a la hora de predecir el siguiente día, así para los 10 días. Realmente es como se predice la demanda eléctrica ya que se compra con un día de antelación (posteriormente hay ajustes). Sin embargo, por falta de tiempo no se ha podido estudiar bien esta alternativa.